#!/usr/bin/perl
# 2012/12/06 Count total observed number of each fragments to define un-alingnable fragments

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 2 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file]\n";
}

my %opt;
getopts("i:", \%opt);
my $FILE_map = $opt{i};

my $FILE_HINDIII_def = '/cscb/noma_nfs/Project/2012-11-15_HiD_1%4%/out/HindIII_sites.txt';


### read Hind III sites
my %name2num;
my %hinDetail;
{
	my $fh_in = IO::File->new($FILE_HINDIII_def) or die "cannot open $FILE_HINDIII_def: $!";
	while($_ = $fh_in->getline()){
		if(m/^#/){
			next;
		}
		s/\r?\n//;
		my ($number, $chr, $pos, $before, $after) = split /\t/;
		my $numLeft = $number - 1;
		my $numRight = $number + 1;

		# left fragment
		{
			my $nameS = $numLeft . 'R';
			my $nameE = $number . 'L';
			$name2num{$nameS} = $numLeft;
			$name2num{$nameE} = $numLeft;
			my $start = $pos - $before;
			my $end = $pos;
			$hinDetail{$numLeft} = [$chr, $nameS, $start, $nameE, $end, $before];
		}

		# right fragment
		{
			my $nameS = $number . 'R';
			my $nameE = $numRight . 'L';
			$name2num{$nameS} = $number;
			$name2num{$nameE} = $number;
			my $start = $pos;
			my $end = $pos + $after;
			$hinDetail{$number} = [$chr, $nameS, $start, $nameE, $end, $after];
		}

	}
	$fh_in->close();
}



#---------------------------------------
# count number
#---------------------------------------
my %SUM;
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	my $num1 = $name2num{$hinID1};
	my $num2 = $name2num{$hinID2};
	$SUM{$num1}++;
	$SUM{$num2}++;
}
$fh_map->close();


#---------------------------------------
# output
#---------------------------------------
print "chr\thinID1\tstart\thinID2\tend\tlength\tobservedNumber\n";
foreach my $k(sort {$a <=> $b} keys %hinDetail){
	my $value = exists $SUM{$k} ? $SUM{$k} : 0;
	print join("\t", @{$hinDetail{$k}}) . "\t" . $value . "\n";
}







